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ABSTRACT 

We carry out a linear analysis of the vertical normal modes of axisymmetric pertur- 
bations in stratified, compressible, self-gravitating gaseous discs in the shearing box 
approximation. An unperturbed disc has a polytropic vertical structure that allows 
us to study specific dynamics for subadiabatic, adiabatic and superadiabatic vertical 
stratifications, by simply varying the polytropic index. In the absence of self-gravity, 
four well-known principal modes can be identified in a stratified disc: acoustic p-modes, 
surface gravity f-modes, buoyancy g-modes and inertial r-modes. After classifying and 
characterizing modes in the non-self-gravitating case, we include self-gravity in the 
perturbation equations and in the equilibrium and investigate how it modifies the 
properties of these four modes. We find that self-gravity, to a certain degree, reduces 
their frequencies and changes the structure of the dispersion curves and eigenfunc- 
tions at radial wavelengths comparable to the disc height. Its infiuence on the basic 
branch of the r-mode, in the case of subadiabatic and adiabatic stratifications, and 
on the basic branch of the g-mode, in the case of superadiabatic stratification (which 
in addition exhibits convective instability), does appear to be strongest. Reducing the 
three-dimensional Toomre's parameter Qzd results in the latter modes becoming un- 
stable due to self-gravity, so that they determine the onset criterion and nature of the 
gravitational instability of a vertically stratified disc. By contrast, the p-, f- and con- 
vectively stable g-modes, although their corresponding uj"^ are reduced by self-gravity, 
never become unstable however small the value of Q^d- This is a consequence of the 
three-dimensionality of the disc. The eigenfunctions corresponding to the gravitation- 
ally unstable modes are intrinsically three-dimensional. We also contrast the more 
exact instability criterion based on our three-dimensional model with that of density 
waves in two-dimensional (razor-thin) discs. Based on these findings, we comment on 
the origin of surface distortions seen in numerical simulations of self-gravitating discs. 

Key words: accretion, accretion discs - gravitation - hydrodynamics - instabilities 
- convection - (stars:)planetary systems: protoplanetary discs - turbulence 



1 INTRODUCTION 

Self-gravity plays an important role in a variety of as- 
trophysical systems. It is a main agent determining the 
dynamical evolution of star clusters, galaxies, various 
types of accretion discs, etc. Particularly in protoplanetary 
discs, that are the central subject of our study, self-gravity 
provides one of the main source of outward angular mo- 
mentum transport thro u gh the e xcitation of density waves 
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Linear stability analysis in a vast majority of cases is 
restricted, for simplicity, to razor-thin, or two-dimensional 
(2D) discs that are obtained by vertically averaging all quan- 
tities. In other words, perturbations are assumed to have 
large horizontal scales compared with the disc thickness. In 
this limit, the well-known Toomre's parameter 

c.r2 
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controls the stability of self-gravitating discs l|Toomrelll963 ). 
In this density wave, which is the only mode in 

a 2D disc, is influenced by self-gravity and thus can be- 
come unstable, as the local dispersion relation for the 
latter clearly dem onstrate s ifColdreich fc Tremaind 1 19781 : 
iBinnev fc Tremaine 1987 ; Bertin et al.lll989t ). 

Stability analysis in a more realistic case of self- 
gravitating three-dimensional (3D) discs is more compli- 
cated. The disc is vertically stratified due to both its own 
self-gravity and the vertical component of the gravity of a 
central object. Depending on the nature of the stratification, 
there exists a whole new set of various vertical modes in the 
system (see below) , some of which can become unstable due 
to self-gravity on horizontal length scales comparable to the 
disc thickness. In this situation, the vertical variation of per- 
turbations is important and for a correct characterization of 
the gravitational instability it is necessary to introduce an- 
other parameter not involving height- dependent variables, 
such as the sound speed in Toomre's parameter. Further- 
more, not all types of stratification permit two-dimensional 
modes, that is, modes with no vertical motions commonly 
occurring in the 2D treatment. For example, in non-self- 
gravitating discs with polytropic vertical stru cture, there 
are no 2D modes jLin et al.lll99(]| : iLubow fc Pri naic 1993E; 
iKorvcanskv fc Pringle 199E ) implying that the dynamics 
does not always reduce to that of the 2D case. Therefore, 
a more accurate stability analysis of self-gravitating discs 
should necessarily be three-dimensional. 

Obviously, before studying the gravitational instability 
of stratified discs, one must first classify and characterize 
vertical normal modes of perturbations in the simplified case 
of no self-gravity. Analysis of the modal structure of strati- 
fied, polytropic, compressible , non- self- grav i t at ing discs has 
been done in several papers: i Ruden et al.l (| 19881 . hereafter 
RPL) IKorvcanskv fc Pringld (|l995l . hereafter KP). lOgilvid 
l|l998h . In convectively stable discs, i.e., with subadiabatic 
vertical stratification, four principal types of vertical modes 
can be distinguished. These modes are: acoustic p-modes, 
surface gravity f-modes, buoyancy g-modes and inertial 
r-modes. The modes are named after their correspond- 
ing restoring forces, which can be well identified for each 
mode at horizontal wavelengths smaller than the disc height 
and are provided by one of the following: compressibil- 
ity/pressure, displacements of free surfaces of a disc, buoy- 
ancy due to vertical stratification and inertial forces due to 
disc rotation, respectively, for the p-, f-, g- and r-modes. In 
the case of superadiabatic stratification, the r- and g-modes 
merge and appear as a single mode, which becomes con- 
vectively unstable for horizontal wavenumbers larger than 
a certain value (RPL); the p- and f-modes remain qualita- 
tively unchanged. For neutral/adiabatic stratification, buoy- 



ancy is absent and the g-mode disappears. The main pur- 
pose of this paper is to investigate how self-gravity modifies 
the frequencies and the structure of the eigenf unctions of 
these modes, which mode acquires the largest positive growth 
rate due to self-gravity and, therefore, determines the on- 
set criterion and nature of the gravitational instability of a 
stratified disc. So, the mode dynamics in the 3D case can ap- 
pear more complex than that in the 2D one, where only the 
density wave m ode can be subject to grav i tationa l instabil- 
ity. Previously, iGoldreich fc Lvnden-Belll (|l965al . hereafter 
GLB) considered gravitational instability in a uniformly ro- 
tating gaseous slab with an adiabatic vertical stratification, 
thereby leaving out all modes associated with buoyancy. 
Other studies also considered the gravitational instability 
of 3D galactic discs, however, the analysis was essentially 
2D, finite-thickness effects were only taken into account by 
mean s of various reduction factors in 2D dispersion rela- 
tions (|Shulll968l : iRomeoll 19921 . Il994l ) . In all these studies, as 
in GLB, the main focus was on finding the criterion for the 
onset of gravitational instability, so that a full analysis of 
various types of vertical normal modes existing in stratified 
self-gravitating discs was not carried out. Actually, we gen- 
eralize the study of GLB to subadiabatic and superadiabatic 
stratifications having different modal structure. 

Another motivation for our study is that the f-mode 
is thought to play an important dynamical role in self- 
gravitating discs. The non-linear behaviour of 3D perturba- 
tions involving large surface distortions, as seen in numer- 
ical simulati ons, has been attri buted to the surface grav- 
ity f-mode l|Pickett et all I2OO0I '). However, this was done 
without analysing the behaviour of other vertical modes 
under self-gravity. It was shown that the f-mode leads to 
a large energy dissipation in the vicinity of the disc sur- 
face, which may facilitate disc cooling, because the energy 
is deposited at smaller optical d epth where it can bo radi- 
ated away rnore q uickly (see e.g.. [johnson fc Gammidl2003l : 
iBolev et al.]|2006l ). Later it was realized that in fact the non- 
linear vertical motions in self-gravitating discs can be much 
more complex than ju st the f-mo de and can have a shock 
character (shock bores. iBolev fc D uriscn 2006). Thus, in the 
3D case, the dynamics of self-gravitating discs is much richer 
and diverse than that of idealized 2D ones and requires fur- 
ther study. To fully understand the origin of such three- 
dimensional effects and what type of instability they are as- 
sociated with, one should start with a rigorous linear study 
of the characteristic properties of all the types of vertical 
normal modes mentioned above, not only the f-mode, in the 
presence of self-gravity. The present work is just a first step 
in this direction. 

Numerical simulations of self-gravi tating discs are of- 
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ten i n the context of g l obal d i scs (e.g., Pickett et al.l 
20031: iRice et al' '2003, '20051: iLodato fc Ricd |2004 
Bolev et al . 2006; Bolcy 200^ and, therefore, are not al- 



20051 : 



ways able to well resolve vertical motions, which, as shown 
in the present study, inevitably arise during the develop- 
ment of the gravitational instability associated with intrin- 
sically three-dimensional modes. So, these simulations may 
not quite accurately capture all the subtleties of the gravita- 
tional instability in 3 D discs. I n this connection, we should 
mention the work by iNelsorJ l|2006l ) that extensively dis- 
cusses the issue of vertical resolution and its importance 
in the outcome of the gravitational instability in numerical 
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simulations of self-gravitation discs. Resolving and analysing 
vertical motions are also crucial for properly understanding 
cooling processes in discs and, particularly, whether convec- 
tion is able to provide sufficiently efi'ective cooling for disc 
fragmentation to occur, wh i ch is still a matter for debate 
in the literature ||BossI|2004| : iMaver et al]|2007l : iBolev et af] 
l2006l . l2007l : [Rafikovl2007f ). In addition, these studies, for sim- 
plicity, use the criterion for gravitational instability based 
on the two-dimensional Toomre's parameter, which, as we 
will demonstrate, cannot always be uniquely mapped into an 
analogous three-dimensional stability parameter and give a 
precise criterion for the onset of gravitational instability. 
In this paper, following other works in a simi- 



lar vein: iLubow fc Pringlel l|l993bl . hereafter LP), KP, 



iLubow fc Qgilviel l|l998l ."hereafter L098) 



we adopt the 

shearing box approximation and consider the linear dynam- 
ics of vertical normal modes of perturbations in a compress- 
ible, stratified, self-gravitating gaseous disc with Keplerian 
rotation. In the unperturbed disc, pressure and density are 
related by a polytropic law, which is a reasonably good ap- 
proximation for optically thick discs (see e.g., L098). This 
allows us to consider the specific features of the dynam- 
ics for subadiabatic, adiabatic and superadiabatic vertical 
stratifications by simply varying the polytropic index. As a 
first step towards understanding the effects of self-gravity 
on the perturbation modes, we restrict ourselves to axisym- 
metric perturbations only. The linear results obtained here 
will form the basis for studying the non-linear development 
of gravitational instability in the local shearing box approx- 
imation that allows much higher numerical resolution than 
global disc models can afford. 

The paper is organized as follows. The physical model 
and basic equations are introduced in Section 2. The clas- 
sification of vertical modes in the absence of self-gravity is 
performed in Section 3. Effects of self-gravity on all normal 
modes are analysed in Section 4. In Section 5, we focus on 
the properties of gravitational instability in 3D. Compari- 
son between the criteria of gravitational instability in 3D 
and 2D is made in Section 6. Summary and discussions are 
given in Section 7. 



2 PHYSICAL MODEL AND EQUATIONS 

In order to study the dynamics of three-dimensional 
modes in gaseous self-gravitating discs, following LP, 
KP, L098, we adopt a local shea ring box approxima- 
tion iGoldreich fc Lvnden-Be"illl965bl '). In the shearing box 
model, the disc dynamics is studied in a local Cartesian ref- 
erence frame rotating with the angular velocity of disc ro- 
tation at some fiducial radius from the central star, so that 
curvature effects due to cylindrical geometry of the disc are 
ignored. In this coordinate frame, the unperturbed differ- 
ential rotation of the disc manifests itself as a parallel az- 
imuthal flow with a constant velocity shear in the radial di- 
rection. A Coriolis force is included to take into account the 
effects of coordinate frame rotation. The vertical component 
of the gravity force of the central object is also present. As 
a result, we can write down the three-dimensional shearing 
box equations 



du 

H 



- 20z X u + V 



(^-Qf7V + ^n^z' + -Vp = 0, (1) 



and the equation of conservation of entropy 
dt {p-f J 



(2) 



(3) 



This set of equations is supplemented with Poisson's equa- 
tion to take care of disc self-gravity 



+ 



dx^ dy' 



+ 



(4) 



Here u — (u^ ,Uy,Uz) is the velocity in the local frame, p, Pjip 
are, respectively, the pressure, density and the gravitational 
potential of the disc gas. is the angular velocity of the local 
reference frame rotation as a whole, x, y, z are, respectively, 
the radial, azimuthal and vertical coordinates, z is the unit 
vector along the vertical direction and d/dt = d/dt + u -V 
is the total time derivative. The shear parameter is g = 
1.5 for the Keplerian rotation considered in this paper. The 
adiabatic index, or the ratio of specific heats, 7 = 1.4 as 
typical of a disc composed of molecular hydrogen; we adopt 
this value throughout the paper. 



2.1 The equilibrium disc model 

Equations (1-4) have an equilibrium solution that is station- 
ary and axisymmetric. In this unperturbed state, the veloc- 
ity field represents, as noted above, a parallel azimuthal fiow, 
uq, with a constant radial shear q: 

UxO = ItzO = 0, UyO = ~qQx. 

In the shearing box, the equilibrium density po, pressure po 
and gravitational potential ipo depend only on the vertical 
z— coordinate and satisfy the hydrostatic relation 

1 dpo 
po dz 



90 = 



n z + — , 
dz 



(5) 



dz^ 



= AnGpo- 



(6) 



As in LP, KP and RPL, pressure and density in the unper- 
turbed disc are related by a polytropic relationship of the 
form 



Po 



(7) 



where K is the polytropic constant and s > is the poly- 
tropic index. The Brunt- Vaisala frequency squared is defined 



po \ci dz dz J \ 



s + 1 



II 50 

ci' 



(8) 



where = jpo/po is the adiabatic sound speed squared. 
If 1 -f 1/s < 7 (subadiabatic thermal stratification), then 
Nq > all along the height and the equilibrium vertical 
structure of the disc is convectively stable. In the oppo- 
site case 1 -I- 1/s > 7 (superadiabatic thermal stratifica- 
tion) , A'^ < everywhere and this corresponds to a convec- 
tively unstable equilibrium. For 14- 1/s = 7 (adiabatic ther- 
mal stratification), A^q = and all the motions/modes due 
to buoyancy disappear. Later we will consider these three 
regimes separately. 
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To determine the vertical structure, we need to solve 
equations (5-6) subject to the boundary condition that the 
pressure vanish at the free surface of the disc. Because we 
have a polytropic model, it is convenient to work with the 
pseudo-enthalpy 



wo 



l/s 



The disc structure is also symmetric with respect to the mid- 
plane, 2 = 0, and, as a consequence, it follows from equa- 
tions (5-7) that the derivative of wo at the midplane vanishes 
(i.e., dwo/dz{0) — 0). Because of this reflection symmetry, 
we consider only the upper half of the disc, z J5 0. At the 
surface of the disc wq = 0, similar to the pressure. This al- 
lows us to determine the equilibrium vertical structure of the 
disc. The non-dimensional variables being used throughout 
the paper are: 



Po 



Csm 



y 



wo 





zQ. 

Z -> , 






Wo 









where = po(0), Wm = ■!i'o(0), Csm = £^(0) are the mid- 
plane values of the equilibrium density, pseudo-enthalpy and 
sound speed. We define the three-dimensional analogue of 
Toomre's parameter as 



47rGp„ 



which is a measure of disc self-gravity (from now on until 
Section 6, we will use Qzd without subscript everywhere, so 
it should not be confused with the 2D Toomre's parameter). 
Using that at the midplane dwo/d2(0) = 0, we can derive 
from equations (5-7) a single equation for the normalized 
pseudo-enthalpy 



s + 1 ldwo\ 
27 \ dz j 



(1 - wo) 



Q(s+ 1) 



(i-^r^). 



(9) 



The normalized density and sound speed are found from 
Po = Wq, Cs = Wo- Equation (9) shows that even though 
the disc is in Keplerian rotation (i.e., self-gravity in the ra- 
dial direction can be neglected), it must be included in de- 
termining the vertical structure. We integrate equation (9), 
starting at z = with «;o(0) = 1, until wo, monotonically 
decreasing with 2;, reaches zero at some finite height/edge 
z = h, which is therefore determined as a result of the inte- 
gration process. At this edge, the density and sound speed 
also vanish, po{h) — 0,Cs{h) = (see also GLB and RPL). 
The height and entire vertical structure of a polytropic disc 
are therefore uniquely specified solely by Q and s, which are 
free parameters in equation (9). 

Figure 1 illustrates the sub- and superadiabatic equi- 
librium vertical structures obtained from equation (9) for 
various Q. For a fixed s, the disc height h decreases with 
decreasing Q. The sound speed is a decreasing function of z 
that does not permit the existence of two-dimensional modes 
in three-dimensional models of poly tropic discs as opposed 
to isothermal ones (ILin et al. I ll99(il . LP, KP). We also see 
that the equilibrium structures for sub- and superadiabatic 
cases look similar except that they have, respectively, every- 
where positive and negative Ng that diverges at the surface 
because the sound speed vanishes there. 






Figure 1. Vertical variations of the normalized density and 
squared sound speed for subadiabatic {s = 5) and superadia- 
batic (s = 1.7) equilibrium states with Q = 0.1 (thick solid lines), 
Q = 0.5 (dashed-dotted lines), Q = 1 (dotted lines), Q = 10 (thin 
solid lines) and the non-self-gravitating case {Q 00, dashed 
lines). The lines corresponding to Q = 10 are close to those in the 
non-self-gravitating limit, implying that the role of self-gravity is 
negligible for Q ^ 10. The disc heights, h, are different depending 
on Q and s. 



2.2 Perturbation equations 

We consider here small axisymmetric {d/dy = 0) perturba- 
tions of the form u', p',p', 1/)' cx exp[—iujt + ikx) about the 
equilibrium states found above, where ui and k are the fre- 
quency and the radial wavenumber, respectively. Without 
a loss of generality, we assume throughout that k is non- 
negative, 0. After switching to non-dimensional vari- 
ables: 



kCsm J. J No 



50 



I 

p 

pm ' 



y 



P 



2 ' ^ 2 ' 

PmCsm Csm 



we can derive from equations (1-4) the following set of equa- 
tions governing the linear dynamics of axisymmetric pertur- 
bations 



■ ' , o ' -11' 

luju^ — h 2ij„ — ikih 

po 

iuju'y = {q~ 2)u'^ 

. , 1 dp' p' dtp' 

iu)u-^ = ; go — 

po dz po dz 

- iujp + ikpou'^ + '^(Po^'z) = 



iuj{p' - CsP ) + 



2 , PoCsNo , 



go 



= 



(10) 

(11) 

(12) 
(13) 
(14) 
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Q' 



(15) 



p = goUz at z = h. 



(21) 



If we now make the changes: — >■ u^, up' /po — >■ p', lj?/;' — >■ 
•0' and ehminate Uy, p' from equations (10), (11), (13), 
we arrive at the following set of equations for the three basic 
quantities u'^, p', ■(/)' (henceforth primes will be omitted) 



dz 

dp 
dz 



90. 



P + [ijJ 

90 



P + 



dz 



dz'2 



Q 



P ^ 
ci go 



(16) 
(17) 
(18) 



where the non-dimensional epicyclic frequency k is given 
by = 2(2 — q) = 1. Notice that the density pertur- 
bations on the right hand side of the linearized Poisson's 
equation (18) consist of two physically different parts. The 
density perturbations due to compressibility, which are in- 
versely proportional to the squared sound speed and pro- 
portional to the pressure perturbation, and the density per- 
turbations due to the stratified background, which are pro- 
portional to the vertical velocity perturbation and A^o • This 
latter term due to stratification is obviously absent in the 
two-dimensional analysis of gravitational instability as well 
as when the disc is adiabatic (see GLB). Equations (16- 
18) form the basis for determining the axisymmetric normal 
modes of perturbations existing in compressible, stratified, 
polytropic, self-gravitating discs. In the non-self-gravitating 
limit, these equations reduce to the main equations of KP. 
The reflection symmetry of the unperturbed disc allows us 
to take the eigenfunctions/modes to be either even or odd 
with respect to z. So, to determine the eigenfrequencies (dis- 
persion diagrams, that is, uj as a function of k for various 
branches) and eigenfunctions, we can numerically integrate 
the main equations (16-18) only in the upper half, !^ z ^ h, 
of the full vertical extent of the disc provided that suitable 
boundary conditions are imposed at 2; = and z = h. In 
other words, we need to solve a boundary value problem. 
Here the disc height h, as noted above, is determined by 
the parameters Q and s. By inspecting equations (16-18), 
it is clear that a normal mode whose pressure p and po- 
tential V perturbations are even functions under reflection 
in 2, has odd vertical velocity and gravitational poten- 
tial derivative dtp/dz perturbations under reflection and vice 
versa. We define a normal mode as being 'even' ('odd') if 
the corresponding vertical velocity and gravitational poten- 
tial derivative are o dd (even) func tions of z (this convention 
agrees with that of IOgilvi3 Il998l . but differs from LP and 
KP). Thus, the boundary conditions at the midplane, z = 0, 
for the even modes are: 



u40) = 0, di>/dz{0) = 
and for the odd modes are: 
p(0) = 0, ^{0) = 0. 



(19) 



(20) 



At 2 = /i we impose the usual free-surface boundary con- 
dition for which the Lagrangian pressure perturbation van- 
ishes. In our new non-dimensional variables this translates 



The boundary condition for the gravitational potential can 
be derived by treating the surface displacement as a surface 
distribution of gravitating matter at 2 = /i. Using the conti- 
nuity of potential across the boundary. Gauss's flux theorem 
and the fact that outside the disc, at 2 — >■ ±00, gravitational 
potential should vanish, we arrive at the following condition 
at the disc surface (see e.g., GLB. lLubow fc Pringle|[l"993al ) 



dz 



+ kip ■ 



Q 



at 



h. 



(22) 



Equations (16-18) supplemented with boundary conditions 
(21), (22) at the surface 2 = /i and (19), (20) at the mid- 
plane fully determine a boundary value problem. Using these 
boundary conditions and variational principle, it can be eas- 
ily shown that lj^ is real (GLB, RPL) that much alleviates 
the search of eigenfrequencies. We integrate these equations 
from z = h downwards to 2 = separately for the even 
modes under conditions (19) and for the odd modes un- 
der conditions (20). We will see below that in the presence 
of self-gravity, the dispersion diagrams for these two mode 
parities are quite different. In other words, self-gravity in- 
fluences even and odd modes, or changes their dispersion 
characteristics, in different ways. As is typical, for a given 
equilibrium structure, that is, for given Q and s, and for a 
given radial wavenumber k, midplane boundary conditions 
can be satisfied only for certain values of uj, yielding disper- 
sion relations uj{k) for different branches of modes classified 
below. We use a Runge-Kutta integrator and root-finding 
algorithms of MATLAB package to first numerically solve 
differential equations (16-18) and then find eigenfrequencies. 



3 THE CLASSIFICATION OF VERTICAL 
MODES 

In this section, we for the moment remove self-gravity (put 
V' — > 0, Q — >■ 00) from the perturbation equations (16-18) 
and from the equilibrium in order to classify and character- 
ize all the vertical axisymmetric normal modes existing in a 
three-dimensional disc. In special cases, an analogous classi- 
fication of modes in non-self-gravitating discs has been don e 
previously by several authors (RPL, LP, KP. IOgilvij|l998l ). 
The aim here is to briefly review and synthesize the results 
from these studies. In the next section, we will again switch 
on self-gravity and investigate how it affects the characteris- 
tics of these modes. Accordingly, in the non-self-gravitating 
case, we need the midplane and surface boundary conditions 
discussed above, but only for pressure and vertical velocity 
perturbations. 

Figure 2 shows the typical dispersion diagrams for three 
types of equilibria: subadiabatic (1 -I- 1/s < 7), adiabatic 
(1 + 1/s = 7) and superadiabatic (1 + 1/s > 7) stratifica- 
tions. In the subadiabatic case, one can identi fy four distinc t 
classes of vertical normal modes (see also KP, IOgilvieiri998l l . 
These modes are: the acoustic p-modes, the surface grav- 
ity f-modes, the buoyancy g-modes and the inertial r-modes 
the restoring forces for which at large radial wavenumbers 
kh ^ 1 are mainly provided, respectively, by compressibil- 
ity/pressure, by the displacements of free surface of the disc, 
by buoyancy due to the vertical entropy gradient and by 
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Figure 2. Dispersion diagrams in the non-self-gravitating case for a subadiabatic disc with s = 5 (panels a, b, c), for an adiabatic 
disc with s = 2.5 (panels d, e) and for a superadiabatic disc with s = 1.7 (panels f, g). The solid lines in panels (a),(d),(f) show the 
frequencies for different branches of the surface gravity f- and acoustic p-modes vs. radial wavenumber k. In order of increasing frequency 
these branches/curves are f°, f°, pf , pf, p|, p§, P3, P3, p|, pj, P5, P5 (superscripts denote the even and odd modes, which are numbered, 
respectively, according to the number of nodes of vertical velocity and pressure perturbations in the interval < z ^ h). Although, 
for large fc, the frequencies of even and odd modes coalesce. The dashed lines in these panels correspond to frequencies of the f-modes 
computed in the incompressible limit. The solid lines in panel (b) show the convectively stable g-mode branches, which in order of 
decreasing frequency are g°'°, gj''^, gg''^, g^''^, g^'^. Here the frequencies of even and odd modes coincide. The solid lines in panels (c),(e) 
show the r-mode branches, which in order of increasing frequency are rjj, rQ, rj, r°, rj, rj, r^, r^ {r^ and are the basic even and odd 
r-modes and subscript '0' means that, respectively, vertical velocity an d pressure pe rturbations for them have no nodes in the interval 
< z ^ /i. In this respect, our numbering of modes differs from that of IOeilvi3ll99^ . The solid lines in panel (g) show the convectively 
unstable g- modes, which in order of increasing cj-^ are g°, gp, gj, g|, gg, gj, g4, gg (gg denotes the basic even g-mode, the corresponding 
vertical velocity perturbation for which has no nodes in the range < z ^ h). Similar to the p-modes, the convectively unstable even and 
odd g-modes coalesce at large k. In each panel, the dashed lines show the corresponding mode branches computed for the incompressible 
case with the same ordering of eigenfrequencies. From panels (b),(g), we see that compressibility most strongly affects the convectively 
stable and unstable g-modes. 



inertial forces due to disc rotation. The existence of the g- 
modes in polytropic discs is attributable to the variation of 
the sound speed and No with height; the latter diverges at 
the disc surface giving rise to the trapped g-mode there. For 
1 4- 1/s < 7, all these modes have uj^ > and, therefore, 
are stable. The p-, f- and g-modes are high-frequency modes 



with frequencies always larger than the epicyclic frequency, 
uj^ K^, while the r-mode is of low- frequency with cj^ ^ 
(we remind that hereafter k = 1). 

For adiabatic stratification, there is no buoyancy (A'^^ — 
0) and, therefore the g-mode disappears, while other modes 
remain qualitatively unchanged. Similarly, the r-mode does 
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not exist in a disc that has zero epicyclic frequency. For 1 + 
1/s > 7, the g-mode becomes convectively unstable (u^ < 0) 
for radial wavenumbers larger than a certain value because 
of negative buoyancy. Its characteristic timescale (growth 
rate) can be of the order of the epicyclic frequency or less 
and due to this, it interferes with motions corresponding to 
the r-mode. Consequently, in the superadiabatic case, the r- 
and g-modes merge at uj^ ^ and appear as a single mode, 
which we still call the g-mode and not the r-mode, because 
the behaviour of corresponding eigenfunctions with height in 
this case is similar to that of the convectively stable g-mode 
(see also RPL); the p- and f- modes are not much affected. 

All these modes come in even and odd pairs. In Fig. 2, 
even (odd) modes are numbered according to the number 
of nodes of the vertical velocity (pressure) perturbation in 
the interval < z ^ h (the node at 2; = is not counted). 
We also see the coalescence of the dispersion curves of the 
even and odd p- and f-modes and also the convectively un- 
stable even and odd g-modes as k increases past a certain 
point, which depends on the mode number (RPL, KP). This 
is associated with a transition from oscillatory to evanes- 
cent behaviour of eigenfunctions. We do not plot eigenfunc- 
tions here, as an extensive discussion of thei r prop erties 
can be found in RPL, KP, L098 and lOgilvid (|l998t ). We 
only mention that the eigenfunctions of the p- and g-modes 
are trapped near the surfaces and decay towards the mid- 
plane, while the eigenfunctions of the r-mode are concen- 
trated near the midplane and decay towards the surfaces. 
As for the eigenfunctions of the fundamental f-mode, they 
have no nodes and monotonically decay from the surfaces 
to the midplane. Below we show that these properties of 
eigenfunctions are altered due to self-gravity. Specifically, 
the number of nodes along each branch of the dispersion di- 
agrams, which is preserved in the non-self-gravitating limit, 
is no longer constant in the presence of self-gravity. The fre- 
quencies of the p- and r-modes increase with mode number. 
The frequencies and growth rates of the convectively sta- 
ble and unstable g-modes, respectively, decrease with the 
mode number. These are the well-known general properties 
of vertical modes in polytropic discs. 

Although each mode type has its dominant restoring 
force, one out of the above mentioned four types (compress- 
ibility, surface gravity, buoyancy, inertial forces), the other 
three forces also contribute to a certain degree for small ra- 
dial wavenumbers kh 1. For example, the p- and f-modes 
are modified both by rotation/inertial forces and buoyancy, 
the f-mode is also modified by compressibility, the g-mode 
is modified by rotation and compressibility and the r-mode 
is modified by buoyancy and compressibility. For the sake 
of comparison with the relevant result from previous stud- 
ies, we are particularly interested in to what degree com- 
pressibility modifies generally non-compressive f-, g- and r- 
modes. So, we decided to juxtapose the dispersion diagrams 
for these modes computed separately in the compressible 
and incompressible cases. We take the incompressible limit 
by formally let ting the adiabatic index go to infinity, 7 — >■ 00 
l|Ogilvielll99i '). After that equations (16-18) without self- 
gravity take the form 

duz k^ 



We solved these equations with the same boundary condi- 
tions and found the dispersion diagrams of the f-, g- and 
r-modes that survive in this limit (obviously, the acoustic 
p-mode disappears). The results are plotted in Fig. 2 with 
dashed lines. Notice that by taking the incompressible limit 
in this manner, we have been able to retain the f-mode, as ex- 
pected. Instead, KP set the density perturbations to zero in 
the linearized continuity equation (anelastic approximation) 
that resulted in the f-mode disappearing in their incompress- 
ible limit. We see that compressibility most strongly affects 
the convectively stable and unstable g-mode branches with 
small mode numbers (Figs. 2b, 2g), or equivalently with ver- 
tical extent compara ble to the di sc height, even at large 
radial wavenumbers (|Ogilvielll99"8h . The frequencies of the 
f- and r-modes do not change much, indicating that these 
modes are nearly incompressible. 



4 VERTICAL MODES IN THE PRESENCE OF 
SELF-GRAVITY 

Having characterized all types of axisymmetric normal 
modes in a stratified disc, let us now compute the disper- 
sion diagrams taking into account self-gravity in the per- 
turbation equations and in the equilibrium. This will allow 
us to understand how the frequencies and the structure of 
the eigenfunctions of the above-described mode types are 
altered by self-gravity, which mode acquires the largest pos- 
itive growth rate in the presence of self-gravity and, thus 
determines the onset criterion and nature of gravitational 
instability of a disc. In other words, we return to the bound- 
ary value problem formulated in Section 2, which is repre- 
sented by equations (16-18) supplemented with boundary 
conditions (21), (22) at the surface and conditions (19), (20) 
of the even and odd symmetry of a solution at the midplane. 
For each Q and s, we first determine the corresponding ver- 
tical distribution of the equilibrium quantities in equations 
(16-18), po, c^, go, Nq , with height as described in Section 2 
and then based on this compute the normal modes. 

Figures 3,4,5 and 6 show the typical dispersion dia- 
grams for the p-, f-, g- and r-modes in a self-gravitating 
disc for subadiabatic, adiabatic and superadiabatic vertical 
stratifications. We separately plot the dispersion diagrams 
of the even and odd parities for each mode type to clearly 
see the influence of self-gravity on them, which, as evident 
from these figures, depends on the mode parity. Unlike in 
the non-self-gravitating limit, the number of nodes of the 
vertical velocity and pressure perturbations in the presence 
of self-gravity are not preserved along each mode branch. 
However, as we will see below, at large k, the infiuence of 
self-gravity on the mode dynamics is small and the disper- 
sion diagrams merge with their non-self-gravitating coun- 
terparts (shown with dashed curves in Figs. 3, 5). So, the 
naming of the modes for smaller k, where the effect of self- 
gravity is important, is done by continuity with the large-fe 
limit for each mode branch. 
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Figure 3. Dispersion diagrams of vertical modes in a subadiabatic self-gravitating disc with s = 5 and Q = 0.2. Shown are (a) the even 
p- and f-modes, (b) the even g-mode, (c) the even r-modc, (d) the odd p- and f-modes, (e) the odd g-modc, (f) the odd r-mode. Dashed 
lines are the corresponding branches computed without self-gravity in the perturbations, but with the same self-gravitating equilibrium. 
The frequency ordering and mode naming are the same as in Fig. 2. Self-gravity reduces/shifts the frequencies of the even and odd p-, 
f-modes and the even g- and r-modes mainly in the range ^ A; ^ 2; the frequencies of the odd g- and r-modes are almost unaffected by 
self-gravity. The frequency of the basic even r-mode (denoted above as rg) is modified most strongly by self-gravity as evidenced by the 
largest dip on the corresponding dispersion curve in panel (c). For large k and/or large mode numbers, the effect of self-gravity becomes 
weaJj and the dispersion curves merge with the dashed ones for the non-self-gravitating case. 



4.1 Subadiabatic stratification 

Consider first the subadiabatic case (Fig. 3), where we again 
have the p-, f-, g- and r-modes, but their dispersion diagrams 
are modified/shifted from their non-self-gravitating counter- 
parts towards lower values due to self-gravity. As illustrated 
in Figs. 3a, 3d, self-gravity reduces the frequencies of all 
branches of the p- and f-modcs, for both oven and odd pari- 
ties, but it more affects the f-rnodos. The situation for the g- 
and r-modes is different (Figs. 3b,3c,3e,3f): only the frequen- 
cies of the even g- and r-modes are reduced by self-gravity 
mostly for radial wavcnumbors in the range S5 fc ^ 2 (dips 
on the corresponding dispersion curves in Figs. 3b, 3c indi- 
cating deviations from the non-self-gravitating dashed ones), 
whereas the frequencies of the odd g- and r-modes are al- 
most unaffected by self-gravity. The frequencies of the fun- 
damental f-modcs, the first few branches of the p-niodcs and 
also the first few branches of the even g- and r-modes are 
reduced noticeably. With increasing mode number amd/or 
radial wavenumber, the effect of self-gravity on the eigenfre- 



quencies gradually falls off, because the corresponding eigen- 
functions become of smaller and smaller horizontal and/or 
vertical scale (the number of nodes increases) and take the 
form similar to that in the non-self-gravitating limit. As 
a result, the dispersion diagrams become more and more 
identical to their non-self-gravitating counterparts. The or- 
dering of frequencies for the even and odd p-, f-, g- and 
r-modes, as it is in the non-self-gravitating case (see Fig. 2): 
a;^(r) ^ ^ ^^(g) < '^^(f) < '^^(p)- is preserved in the 
self-gravitating case as well, however small the value of Q. 
Note also in Fig. 4 that the frequencies of the the even and 
odd f-modes and also of the basic even g-mode, for which 
the influence of self-gravity is stronger, never fall below the 
epicyclic frequency with decreasing Q and, therefore, these 
modes always remain gravitationally stable (we also checked 
that the same situation holds for the dispersion curves of the 
even and odd f-modes in the adiabatic and superadiabatic 
cases described below). This is a consequence of the three- 
dimensionality of the disc. Thus, as evident from Figs. 3c, 4d, 
self-gravity most significantly affects the basic branch of the 
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4.2 Adiabatic stratification 

In the adiabatic case (Figs. 5a-5d), there is no g-mode and 
the behaviour of the p-, f- and r- modes under self-gravity 
has qualitatively the same character as in the subadiabatic 
case above. Again, the basic branch of the even r-mode ap- 
pears to be most significantly affected by self-gravity (the 
corresponding dispersion curve in Fig. 5c has the same form 
as that in Fig. 3c) and becomes gravitationally unstable. 
The ij^ of the p- and f-modes, although lowered by self- 
gravity, always remain larger or equal to ^ irrespective of 
the value of Q. In this case, s = 2.5 and the corresponding 
Qcr ~ 0.177 and km = 0.81. Again the radial scale of insta- 
bility A,„ = 2n /k,n = 7.75 is not far from the disc thickness 
for these parameters, H — 2h = 2.08. 



Figure 4. Dispersion curves for s = 5 and Q = 0.01 (dashed- 
dottcd lines), Q = 0.05 (dashed lines), Q = 0.2 (solid lines). Pan- 
els (a),(b) show the even and odd f-modes, respectively. Panel (c) 
shows the basic even g-mode and panel (d) shows the basic even 
r-mode. The frequencies of the f- and g-modes never fall below 
the epicyclic frequency with decreasing Q and, therefore, are al- 
ways stable. The basic even r-mode, which is most significantly 
affected by self-gravity, becomes unstable when the dip on the 
corresponding dispersion curve crosses the y = —axis. This dip 
deepens and broadens with decreasing Q. 



even r-mode and only this branch can become unstable due 
to self-gravity. From Fig. 4d, it is seen that the dip on this 
branch deepens and broadens with decreasing Q. (We de- 
note the wavenumber at which the minimum of this dip is 
located by km-) The gravitational instability sets in when 
the dip's minimum first touches the o;^ — 0— axis at some 
Qcr and then farther extends into the unstable oj^ < region 
for Q < Qcr- We will demonstrate in Section 5 that the basic 
even r-mode becomes strongly compressible at fc ~ km (Fig. 
8), as opposed to the r-mode in a non-self-gravitating disc, 
and the density perturbations due to compressibility are re- 
sponsible for its gravitational instability. For s = 5, we find 
Qcr = 0.168, which gives H = 2h = 2.96 for the disc total 
thickness, and km = 0.8 (see Fig. 12). Note that the charac- 
teristic radial scale of instability Xm = 2Ti/km ~ 7.85 is not 
much greater than the disc thickness H - In self-gravitating 
discs, compressibility and inertial forces play an important 
role in the dynamics of the basic even r-mode at wavelengths 
A ~ Xm ~ H and so it resembles a 2D density wave at such 
wavelengths (see Section 4.4). The fact that the character- 
istic scale of gravitational instability Xm turns out to com- 
parable to the disc thickness may offer a clue why angular 
momentum transport due to self-gravity tends to be a lo- 
cal phenomenon (e.g., lGammiell200ll : iLodato fc Ricell2004l : 
iBolev at al.ll200^ ). In this respect, it should also be men- 
tioned that the analogous gravitational instability of low- 
frequency mo des in a strongly compress ed gaseous slab was 
also found bv iLubow fc Pringld l|l993al ). In their analysis, 
these modes, called neutral modes, are basically degenerate 
r-modes, since the compressed gaseous slab was not rotating. 



4.3 Superadiabatic stratification 

The superadiabatic case (Figs. 5e-5h) is interesting, because 
at y ^ y, as classified above, instead of the r-mode there 
is the g-mode, which in addition to convective instability 
can also exhibit gravitational instability. As in the non-self- 
gravitating case, for radial wavenumbers larger than a cer- 
tain value, all the branches of the g-mode become unstable 
because of the negative vertical entropy gradient. From Fig. 
5g we see that for ^ A: 2, self- gravity produces dips on 
the g-mode dispersion curves, again preferably for the even 
parity ones. When Q drops sufficiently, the dip on the ba- 
sic branch of the even g-mode, which appears to be most 
affected by self-gravity, starts to cross the lj^ = 0-axis in a 
way similar to that of the basic even r-mode branch above 
and this signals the onset of gravitational instability. For 
s = 1.7, adopted in Fig. 5, we find Qcr = 0.184, which 
gives a thickness H — 2h = 1.76, and km = 0.83, so that 
Xm ~ 2-K/km = 7.57 H- As for the p- and f-modes, they 
behave under self-gravity in exactly the same manner as in 
the above cases. In particular, although their frequencies 
are reduced, they can never become gravitationally unsta- 
ble even for very small values of Q- Thus, for superadiabatic 
stratification, the basic even g-mode can exhibit two types 
of instabilities: gravitational and convective (see Fig. 6). As 
we see from Figs. 5g,6a, at moderate Q ~ Qcr, the radial 
scales for the activity of self-gravity and convective insta- 
bility are well separated, so that these two effects do not 
interfere with each other in the linear regime. For Q = 0.1 
in Fig. 6a, self-gravity is dominant at ^ k ^ 2.55, while 
convective instability occurs at fc > 4.27, where the effect of 
self-gravity is weak (a similar situation is for Q = 0.2 in Fig. 
5g). However, for very small Q <C Qcr, the radial scales of 
gravitational and convective instabilities overlap (Fig. 6b), 
but the gas motion for such radial scales in the case of very 
strong gravitational instability would hardly resemble con- 
vective motions (it will look more like that in Fig. 11 below); 
convection is simply disrupted by gravitational instability. 
Therefore, we can conclude that unless a disc is strongly self- 
gravitating, self-gravity has little infiuence on the properties 
of convective motions and on the (Schwarzschild) criterion 
for the onset of convective instability that is equivalent to 
iVo^ < 0. 
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Figure 5. Dispersion diagrams for vertical modes in a self-gravitating disc with Q = 0.2 for an adiabatic stratification with s = 2.5 
(panels a-d) and for a superadiabatic stratification with s = 1.7 (panels e-h). Shown are (a,e) the even p- and f-modes, (b,f) the odd p- 
and f-modes, (c) the even r-mode, (d) the odd r-mode, (g) the convectively unstable even g-mode and (h) the odd g-mode. The frequency 
ordering and mode naming are the same as in Fig. 2. As before, dashed curves correspond to the dispersion diagrams calculated without 
self-gravity in the perturbations. The behaviour of the p-, f- and r-modes in the adiabatic case is similar to that in Fig. 3. In the 
superadiabatic case, the dispersion curves of the even g-mode have dips caused by self-gravity in the interval ^ fc 2, while the 
odd g-mode is almost unaffected by self-gravity. The largest dip occurs for the basic branch of the even g-mode (denoted above as gg) 
in panel (g), which appears to be most significantly influenced by self-gravity. For large k and/or large mode numbers, the effect of 
self-gravity becomes weak and the dispersion curves merge with dashed ones for the non-self-gravitating case. Due to the superadiabatic 
stratification, the even and odd g-modes are convectively unstable (i.e., have uj'^ < 0) in the range k > 3.91, where the influence of 
self-gravity is small. 



4.4 Analogy with the 2D density wave theory 



From the three cases considered above, we conclude that 
the basic branch of the low-frequency (cj^ ^ k^) even mode, 
which is the r-mode, in the case of subadiabatic and adia- 
batic stratifications, and the g-mode, in the case of superadi- 
abatic stratification, is most subject to the influence of self- 
gravity. These modes determine the gravitational instability 
in vertically stratified discs. As we have clearly seen in all the 
above cases, the radial scale of the instability is comparable, 
though a bit larger than the disc thickness. Am H, which 
in turn implies that a 2D analysis of the gravitational insta- 
bility is at most marginally valid and the rigorous stability 
(linear) analysis of self-gravitating discs should be three- 
dimensional. However, despite the 2D treatment only being 
marginal, we can still find similarities between our results 



and the 2D density wave theory in thin self-gravitating discs, 
which actually appears to do a decent qualitative job. 

Consider an equivalent zero-thickness disc with the 
sound speed CsmQ and the surface density Eq = podz. 
Then the well-known dispersion relat ion of axisymmetric 2D 
density waves in the thin disc is (e.g.. lGoldreich fc Tremaind 
il97S) 

uj'^ = ~ -p^—k + 1, 

where we have used the same normalization as before and 

^ In fact, there are other options for choosing the sound speed as 
a some kind of height-averaged value. We address this question in 
Section 6, but for the present purpose this does not make much 
difference. 
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Figure 6. Dispersion curves of the basic even g-mode for a su- 
peradiabatic stratification with s = 1.7 for (a) Q = 0.1 and 
(b) Q = 0.01. Dashed lines show the same branches computed 

without self-gravity in the perturbations. For Q = 0.1, the basic 
even g-modc exhibits gravitational instability in the vicinity of 
krn = 1.15 and coiivcctive instability in the range k > 4.27, where 
the influence of self-gravity is weak. For Q = 0.01, the effect of 
self-gravity is much stronger and the radial scales of gravitational 
and convective instabilities overlap in the range 4.25 ^ fc ^ 8.5. 



because of that Q2D = Cam^/nCEo- This dispersion rela- 
tion is a parabola with a minimum at the Jeans wavenum- 
ber kj = I/Q2D. This is the wavcnumbcr, at which the 
effect of self-gravity is most prominent, and if Q2D < 1 
it gives the characteristic scale of gravitational instability. 
At small k <Si kj , 2D density waves are dominated by self- 
gravity and inertial forces, while at large k ^ kj, pres- 
sure/compressibility dominates over self-gravity and density 
waves appear as an acoustic mode. At fc ~ fcj all three fac- 
tors can be important. Let us now look at the dispersion 
curves of the basic branches of the even r- and convectively 
unstable even g-modes in Figs. 3c,5c,5g. They have similar 
parabolic shape in the self-gravity and compressibility dom- 
inated regime at A: ~ fem with the linear phase at smaller 
A; <C fern, where only self-gravity and inertial forces play 
a role. This linear phase at long wavelengths is well repro- 
duced by the 2D dispersion relation. Therefore, we can iden- 
tify the wavcnumbcr km, at which the effect of self-gravity 
on the 3D modes is largest, with the Jeans wavenumber kj. 
In some sense, in the case of instability when Q < Qcr, km 
gives a more accurate value for the radial scale, \m, of the 
gravitationally most unstable mode than that given by the 
Jeans wavenumber k,j in the 2D model. For example, from 
Fig. 12 we find that the vertical structure (po) calculated at 
s = 5 and Qcr = 0.168 gives the corresponding Q213 = 0.76 
and k,j = 1.32, which differs from km = 0.8 found above for 
these parameters. Moreover, we will see in Section 6 that the 
criterion for gravitational instability based on the 3D calcu- 
lations, i.e. Q < Qcr, is more exact than the 2D criterion 
Q2D < 1. 

In the 2D case, as shown above, the Jeans wavcnumbcr 
kj = I/Q2D, implying that it is determined solely by com- 
petition between self-gravity and compressibility/pressure. 
It is now interesting to see how an analogous characteristic 
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Figure 7. Dependence of km on Q for s = 5, which very closely 
follows the power law Q^^/^ (shown with dashed line and scaled 
appropriately) . The gravitational instability sets in at Q = Qcr- = 
0.168 that gives km = 0.8. 



wavenumber km depends on Q in the 3D case. From Fig. 
7 we see that this dependence has a power law character 
Q-'/^ which means that the value of km is again set by 
self-gravity and compressibility, as in the 2D case, so that 
the disc rotation plays a role only in determining lu"^ but not 
km- Indeed, the only lengthscale that can be constructed 
from the sound speed c^m , density pm and gravitational con- 
stant G without the angular velocity, of disc rotation is 
Csm/\/Gpm (analogous to the Jeans length of a collapsing 
3D cloud). If wc normalize it by Csm/f^, we get 2^/tvQ^^^, 
giving the above power law dependence for the correspond- 
ing non-dimensional wavenumber y/TrQ~^^^ to which km is 
proportional. 

Thus, the basic even r-mode or, in the case of superadia- 
batic stratification, the basic even g-mode at k ^ km exhibit 
many of the characteristics of 2D density waves and could be 
regarded as their 3D analogues at such radial wavenumbers 
(see also Section 5). If continued to » km, the density 
wave mode would thus connect up to the large- fc parts of 
these two modes dominated, respectively, by inertial forces 
and by negative buoyancy, instead of merging with the com- 
pressible p-mode, as might seem at first sight from the 2D 
dispersion relation. However, the limit k 3> km, equivalent 
to \ <^ H, means the breakdown of the 2D approximation 
and, therefore, the concept of 2D density waves aX k ^ km 
would not be well-defined. 



5 GRAVITATIONAL INSTABILITY IN 3D: 

PROPERTIES OF THE BASIC BRANCH OF 
THE LOW-FREQUENCY EVEN R-MODE 

In this section wc concentrate only on the properties of the 
gravitationally unstable basic even r-mode. As we have seen, 
the behaviour of the convectively unstable basic even g-mode 
under self-gravity has a qualitatively similar character. 
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Figure 8. Dispersion curves of the basic even r-mode at s = 
5,Q = 0.1 with corresponding h = 1.22. Dotted line corre- 
sponds to the non-self-gravitating case, dashed-dotted - to the 
incompressible case with self-gravity and solid line - to the self- 
gravitating case with compressibility. The influence of self-gravity 
on the basic even r-mode is very weak in the incompressible limit 
and the dashed-dotted dispersion curve almost coincides with the 
dotted one in the non-self-gravitating limit; a slight deviation 
(dip) at small k is due to the stratified background. 



5.1 Effect of compressibility 

As noted above, in the Poisson's equation, the density per- 
turbations giving rise to the gravitational potential pertur- 
bations can be due to compressibility and due to the strat- 
ification of the equilibrium vertical structure. As is evident 
from Fig. 2c, in non-self-gravitating discs, the r-mode is 
nearly incompressible. So, it is interesting to see if it still re- 
mains incompressible in self-gravitating discs and which out 
of these two sources of density perturbations is ultimately 
responsible for the gravitational instability. In order to ex- 
plore this, we again computed the dispersion curve of the 
basic even r-mode in the incompressible limit of equations 
(16-18), which take the form: 



duz 
dz 

dp 
dz 



r(p + V') 



p+ a; 

go 



dip 
dz 



dz'' Q go 



(23) 



(24) 



(25) 



supplemented with the same boundary conditions. Notice 
that on the right hand side of Poisson's equation (25), 
only the density perturbation due to stratification remains. 
In Fig. 8, we compare the dispersion curves of the basic 
even r-mode obtained in the incompressible limit to those 
computed for the compressible self-gravitating and non-self- 
gravitating cases. It is clear that in self-gravitating discs, the 
basic even r-mode becomes compressible at fc ~ km thereby 
providing density perturbation for the gravitational poten- 
tial and, therefore, accounting for the large dip on the dis- 
persion curve. The dispersion curve in the incompressible 
limit deviates only very slightly from that in the non-self- 
gravitating case, because of the second source of density 



Figure 9. Vertical structures of the radial velocity Ux, the ver- 
tical velocity Uz, the pressure p and the density p perturbations 
constituting the eigenfunctions of the gravitationally unstable ba- 
sic even r-mode for s = 5,Q = 0.1, h = 1.22, km = 1 with the 
corresponding largest growth rate '.^min ~ ~0.75. Dashed lines 
show the eigenfunctions (scaled arbitrarily for plotting purposes) 
in the non-self-gravitating case for the same parameters except 
ui-^ = 0.99. Notice that in the non-self-gravitating case, the eigen- 
functions are more concentrated near the midplane, whereas in 
the self-gravitating case they vary over the whole vertical extent. 



perturbation, i.e. stratification, which turns out to be much 
smaller than that associated with compressibility. Such a 
comparison also shows that the primary cause of gravita- 
tional ins tability in incompressibl e discs, as described in 
GLB and lLubow fc Pringld (|l993ah , is the displacements of 
free-surfaces. In our case, the equilibrium density vanishes 
at the surface, so that the effect of surface displacements on 
the growth rate of instability is null (see boundary condition 
22) and in the incompressible case only density perturbation 
due to stratification, as shown, is hardly sufficient for grav- 
itational instability. 



5.2 Eigenfunctions and spatial structure 

Here we compute the vertical structure of the eigenfunctions 
of the gravitationally unstable basic even r-mode and also 
find what type of motions it induces. We take Q — 0.1 and 
s = 5, which gives h = 1.22 for the disc height. For these 
parameters, the dispersion curve of this mode in Fig. 8 has 
a minimum LOmi„ ~ —0.75, i.e. gives the largest growth rate 
of gravitational instability, at km ~ 1- In Fig. 9, we plot the 
corresponding eigenfunctions only in the upper half of disc's 
full vertical extent. Because the mode has even parity, the 
vertical velocity and the derivative of the gravitational po- 
tential are odd functions and, correspondingly, the pressure 
and potential are even functions of z. Only the pressure per- 
turbation has one node in the interval < z ^ h, other quan- 
tities have no nodes and vary over the whole vertical extent. 
This demonstrates that the spatial structure of the gravita- 
tionally most unstable mode is three-dimensional with non- 
zero vertical velocity. To see how self-gravity changes the 
form of the eigenfunctions, we also show the eigenfunctions 
of the same branch in the non-self-gravitating limit, which 
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Figure 10. Eigenfunctions of the even f-mode, which unlike the 
basic even r-mode is gravitationally stable, with the same values 
of parameters as in Fig. 9. The corresponding eigenfrequencies at 
km. = 1 are o;^ = 1.1 for the self-gravitating case (solid lines) and 
cj^ = 1.6 for the non-self-gravitating case (dashed lines). Again, 
non-self-gravitating eigenfunctions have been scaled arbitrarily 
for plotting purposes. 
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Figure 11. Distribution of the density (upper plot) and temper- 
ature (lower plot) perturbations in (x, z) plane constructed from 
the eigenfunctions of the gravitationally unstable basic even r- 
modo in Fig. 9. Superimposed on the density field are the velocity 
vectors of the induced gas flow. 



appear to be more concentrated towards the midplane than 
the self-gravitating ones and have no nodes in the same ver- 
tical range. Since the f-mode plays an important role in the 
angular momentum and energy transport in stratified discs 
(L098), in Fig. 10 we also plot its eigenfunctions in the 
presence of self-gravity for the same parameters. Compar- 
ing Figs. 9 and 10, we clearly see that self-gravity modifies 
the form of the eigenfunctions of both the f- and r-modes. 
The perturbed quantities for both mode types vary with 
height somewhat similarly, especially the vertical velocities, 
which are, however, still at least an order of magnitude less 
than the horizontal ones, so that the motions can be viewed, 
to leading order, as 2D (see also Fig. 11). Thus, the gas mo- 
tion associated with the gravitationally unstable r-mode is 
of similar type to that of the f-mode, which, in turn, implies 
that the r-mode might be as important as the f-mode, i.e., 
might be able to do similar 'dynamical jobs' as the f-mode 
in self-gravitating discs. As mentioned earlier, the non-linear 
behaviour of the 3D perturbations in self-gravitating discs 
has been attributed so lely to the f-mode dynamics without 
analy sing other modes jPickett et al ]|200Gl : lBolev fc DurisenI 
12009 ). 

The perturbed quantities have the form 
f{z)exp{ikmx), f = (ux,Uz,p,p) and from this we can 
construct the spatial picture of the velocity, density and 
temperature perturbations for the gravitationally unstable 
basic even r-mode by taking the real parts. The temperature 
perturbation is found from the ideal gas equation of state 
and is equal to T = — Cgp/po + 7P/po (here the pressure 
perturbation p is as used in the original equations 10-15, 
i.e., before switching to new variables). Figure 11 shows 
the density and temperature fields constructed in this way 
with velocity vectors showing the direction of gas motion 
superimposed on the density field, which in fact resembles 
the classical density profile of a 2D density wave. Near the 
midplane, matter flows (converges), almost parallel to the 



a;— axis, towards high density regions. With increasing z, 
the flow becomes more arc-like, because of the variation of 
the vertical velocity with height. An analogous velocity field 
in the {x,z) plane was also obser ved in s ome non-linear 
simu lations of self-gravitating discs (|Bolev fc Duriscn 200^; 
iBolev et alj|200(^ . which suggests that the latter may be a 
result of the non-linear development of the type of motion 
we see here in the linear regime. We might also speculate 
that in relatively high-mass discs, the non-axisymmetric 
gravitationally unstable basic r-mode will produce similar 
high-density regions, which can be precursors of spiral 
shock fronts (if the disc does not fragment), and arc- like 
motions in the upper layers that can give rise to non-linear 
shock bores inv olving large d istortions of the disc surface 
as described by iBolev fc Durisen (200d). The temperature 
perturbations are fairly non-uniform as well, varying both 
vertically and radially on comparable scales. Obviously, 
larger temperature perturbations correspond to overdense 
regions that are contracting because matter fiows into 
them, and lower - to underdense regions. 



6 STABILITY CRITERIA IN 2D AND 3D 

In the 3D case, the degree and effect of self-gravity is charac- 
terized by Q-iD = / ^TvGpm, where pm is the value of equi- 
librium density at the midplane (see also GLB), which plays 
here a role similar to that of the standa rd 2D Toomre's pa- 
rameter Q2D = c^n/TrGS jToomrdl 19641 ). However, in three- 
dimensional simulations, researchers often still tend to use 
the 2D Toomre's par ameter to characterize the onset of grav- 
itational instability dRice et all l200l l2005l : iLodato fc Ric^ 
l2004l2005l : lBolev et al.ll2006l ). In this section we investigate 
how a more general and exact 3D criterion of gravitational 
instability can be related to the 2D one. In general, these 
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ate as it does not involve uncertainties in how to average the 
sound speed over height. 

In Fig. 12 we also examine the dependence of the crit- 
ical wavenumbers of the gravitational instability on s. We 
see that the 2D Jeans wavenumbers are always larger than 
the actual more exact wavenumber km at the onset of grav- 
itational instability when Q = Qcr,3D- Thus, in the 3D 
case, the critical wavenumbers depend differently on the 
thermodynamics than in the 2D case, where the critical 
kj — I/Q2D = 1 at the onset of instability. 



0.5 L 



Figure 12. Dependence of the critical Qcr,3D and the corre- 
sponding critical Qi cr- 2D a-nd Q2 cr 2D on the polytropic index 
s. Shown also are the k^. at Q^r 3D and the 2D Jeans wavenum- 
bers kji = 1/Qi cr 2D I kj2 = I/Q2 cr 2D as a function of s. 



two parameters are different, because the 2D parameter con- 
tains the sound speed and surface density, which by defini- 
tion do not vary vertically for razor-thin discs, whereas the 
3D parameter does not contain the sound speed and sur- 
face density, but the midplane value of the density. So, if 
we still want to characterize the 3D instability in terms of 
the 2D Toomre's parameter, we need to use some height- 
averaged, or effective sound speed (surface density can be 
uniquely calculated from the vertical den sity distribution) . 
In some sim ulations the midplane values (|Meifa et al.ll2005l : 
I Bole V et al, whi le in others vertically averaged values 

I Rice et al.l 2003l . l2005l ) of the sound speed are used. Here we 
calculate the corresponding critical Qcr.2D for these cases. 

As is well-known, in razor-thin 2D discs, axisymmetric 
perturbat ions (density w aves) are gravitationally unstable if 
Q2D < 1 (|Toomrelll963 ). Let us now find the critical Qcr,3D 
that determines the gravitational instability in 3D, which, 
as found in this study, is associated with the basic branches 
of the even parity low-frequency modes. In Fig. 12, we show 
the critical Qcr.sD as a function of the polytropic index s. 
Given Qcr,3D, we can now find Qi cr,2D defined in terms of 
the midplane sound speed and Q2 cr,2D defined in terms of 
the vertically averaged sound speed: 



Q 



1 cr,2D — 



2Qc 



{2 cr,2D = 



2Qcr,3D Csdz 

h podz 



where po and Cs, as before, are the normalized equilibrium 
density and sound speed the z-dependence of which, as well 
as the value of h, are determined by Qcr,3D and the poly- 
tropic index s. So, the critical Qi cr,2D and Q2 cr,2D calcu- 
lated with the two different methods are different, though 
not far from each other, and less than unity which is their 
critical value in the 2D case. This again confirms the well- 
known result that 3D discs are more stable and to make 
them unstable smaller value of Q2D is necessary when com- 
pared with that required for razor-thin 2D discs (see also 
GLB). But in such cases, using Q3D seems more appropri- 



7 SUMMARY AND DISCUSSIONS 

In this paper, we have analysed the axisymmetric nor- 
mal modes of perturbations in stratified, compressible, self- 
gravitating gaseous discs with subadiabatic, adiabatic and 
superadiabatic vertical stratifications. First, we performed 
a classification of perturbation modes in stratified discs in 
the absence of self-gravity to compare with previous calcu- 
lations. Four well-known main types of modes can be distin- 
guished: acoustic p-modes, surface gravity f-modes, buoy- 
ancy g-modes and inertial r-modes. The restoring forces 
for these modes for large radial wavenumbers are mainly 
provided by one the following: pressure/compressibility, dis- 
placements of the disc surface, buoyancy and inertial forces 
due to disc rotation for the p-,f-,g- and r-modes, respec- 
tively. For smaller wavenumbers, the restoring force for each 
mode is a combination of these forces. In the case of adi- 
abatic stratification, buoyancy is zero and, therefore, the 
g-mode disappears, while other modes remain qualitatively 
unchanged. For superadiabatic stratification, the g-mode be- 
comes convectively unstable and merges with the r-mode, 
so that only a single convectively unstable mode appears in 
the dispersion diagram at lu'^ ^ k^, which we still call the 
g-mode. Due to the refiection symmetry of the equilibrium 
vertical structure with respect to the midplane, each mode 
comes in even and odd pairs. By our terminology, for even 
(odd) modes, pressure and gravitational potential pertur- 
bations are even (odd), while the perturbations of vertical 
velocity and derivative of potential are odd (even) functions 
of the vertical coordinate. After classifying and character- 
izing modes in the absence of self-gravity, we introduced 
self-gravity in the perturbation equations and investigated 
how it modifies the properties of these modes. We found 
that self-gravity, to some extent, reduces the frequencies of 
all normal modes at radial wavelengths comparable to the 
disc height, but its influence on the basic even r-mode, in 
the case of subadiabatic and adiabatic stratifications, and on 
the basic even g-mode, in the case of superadiabatic stratifi- 
cation, appears to be strongest. With decreasing Q3D, these 
modes become unstable due to self-gravity and thus deter- 
mine the gravitational instability of a vertically stratified 
disc. The basic even g-mode also exhibits convective insta- 
bility due to a negative entropy gradient but, unless the disc 
is strongly self-gravitating, these two instabilities grow con- 
currently in the linear regime, because their corresponding 
radial scales are separated. We also obtained the correspond- 
ing criterion for the onset of gravitational instability in 3D, 
which is more exact than the standard instability criterion in 
terms of the 2D Toomre's parameter, Q2D < 1, for axisym- 
metric density waves in razor-thin discs. By contrast, the 
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p-, f- and convectively stable g-modes have their reduced 
by self-gravity, but never become unstable for any value of 
QsD ■ This is a consequence of the three-dimensionality of the 
disc. The eigenfunctions associated with the gravitationally 
unstable modes are intrinsically three-dimensional, that is, 
have non-zero vertical velocity and all perturbed quantities 
vary over the whole vertical extent of the disc. In this regard, 
we would like to mention that resolving the gravitationally 
most unstable mode in numerical simulations thus reduces 
to properly resolving the disc height (together with resolv- 
ing the corresponding radial scale, which, as shown here, 
appears somewhat larg er than the disc height). So, the cri- 
terion of iNelsoiil l|2006h that at least ~ 4 particle smoothing 
lengths should fit into per scale height may apply in three- 
dimensional SPH simulations. This implies a substantially 
larger number of SPH particles per vertical column because 
the disc itself may extend over many scale heights. He also 
shows that a similar criterion applies to grid-based simula- 
tions. 

Here for simplicity, and also to gain the first insight into 
the effects of self-gravity on the vertical modes in stratified 
discs, we considered only axisymmetric perturbations. Non- 
axisymmetric perturbations are dynamically richer, though 
more complicated, because the phenomena induced by Ke- 
plerian shear/differential rotation - strong transient am- 
plification of perturbations and linear coupling of modes 
(not to be confused with non-linear mode-mode interac- 
tions) - come into play for these type of perturbations. Tran- 
sient (swing) amplification of perturbations (density waves) 
has been studied previousl y in razo r- thin 2D approximation 
dGoldreich fc Lvnden-Beiil Il965bl: iGoldreich fc Tremaind 



Il978l : lToomrelll98ll : I Mamatsashvili fc ChageUshvilil |2007| )" 
From the analysis presented here, we may expect that 
in the linear regime, the non-axisymmetric basic even r- 
mode can undergo larger transient amplification due to self- 
gravity than other modes in the disc. This transient am- 
plification of perturbations may be important for explain- 
ing the large burst phases seen in numerical simulations 
at the initial stage s of the d evelopment of gravitational 
instability in discs dRice et al ] [2003. .20051: iLodato fc Ric3 
l2004l2005l : lBolev et al.ll2006l ). So. in this respect one should 
analyse and quantify the transient amplification of non- 
axisymmetric perturbations in stratified 3D self-gravitating 
discs starting with linear theory. As for the linear coupling 
of modes, it was demonstrated that in non-self-gravitating 
stratified discs, Keplerian shear causes rotational (vortex) 
mode perturb ations to couple with and generate g-mode 
perturbations ijTevzadze et al.ll2008l 'l. In the context of 2D 
discs, it was shown that vortex mode perturbations can 
also exc ite density waves d ue to sh ear l|Bodo et aL 20051: 



Mamats ashvili fc Chaeelishv ili 2007; Mamatsashvili fc Ricd 
2009l : iHeinemann fc Papaloizouil2009l ). In the 3D case, there 



are a larger number of modes in a disc and it is quite possible 
that some of them may appear linearly coupled due to shear 
and, therefore, be able to generate each other during evo- 
lution, especially the f-mode and the basic branch of the r- 
mode (because they vary on comparable vertical scales in the 
presence of self-gravity. Figs. 9, 10). Another related prob- 
lem also of interest is the interaction between self-gravity 



and the MRI in magnetized discs (see e.g.. iFromand 120051 ). 
In particular, how self-gravity can modify the growth rates 
of magnetic normal modes responsible for the MRI. Actu- 



ally, this will be the generalization of th e extens i ve ana lysis 
of normal modes in magnetized discs bv lOgilvi3 (|l998l ). 

iLubow fc Qgilvid (11998') showed that in non-self- 
gravitating discs with polytropic vertical stratification, an 
external forcing preferentially excites the f-mode, because it 
has the largest responsiveness to an external driving com- 
pared to other modes. This mode, propagating through a 
disc, results in energy dissipation near the disc surface. 
Based on these results and partly on th e prop erties of f- 
modes in stellar dynamics, iPickett et al.l (|2000l ) identified 
the behaviour of 3D perturbations in self-gravitating discs 
involving large surface distortions and the resulting energy 
dissipation in the upper layers, with f-mode dynamics. The 
self-stimulated potential was thought to play the role of an 
external/tidal force. However, it is not obvious that the ef- 
fect of a self-stimulated potential is the same as that of the 
external potential. In fact, our analysis has revealed that in 
self-gravitating discs, in addition to the f-mode, the r-mode 
can also be dynamically important, because this mode ap- 
pears to be subject to gravitational instability, while the 
f-mode is not. The eigenfunctions of the gravitationally un- 
stable basic even r-mode differ from those of the r-mode in 
non-self-gravitating discs in that they are no longer concen- 
trated near the midplane and behave somewhat similarly to 
the eigenfunctions of the f-mode: they vary over the whole 
vertical extent of the disc and also involve noticeable per- 
turbations of the disc surface. Consequently, like the f-mode, 
the gravitationally unstable r-mode can, in principle, also in- 
duce gas motion causing large surface distortions and resul- 
tant energy dissipation in the upper layers of the disc, which 
is thought to play a role in enhancing disc cooling (because 
the energy is deposited in the upper layers with small optical 
depth, it can be radiated away more quickly and effectively 
cool the disc, but th is is a subject of further study, see e.g., 
I Johnson fc Gammie 2003; Bolcv et al. 2006). Furthermore, 
in the case of non-axisymmetric perturbations, as mentioned 
above, because of shear, the gravitationally unstable r-mode 
can couple with and generate the strong f-mode. So, the 
surface distortions may be caused by a combination of the 
f- and r-modes. In order to explore where dissipation can 
predominantly occur in a se lf-gravitating dis c, one needs 
to generali z e the analysis of iLin et al.l (|l990l ). LP, L098, 
iBate et al.l l|2002l ') on the propagation of waves in stratified 
non-self-gravitating discs and consider the propagation of 
non-axisymmetric modes in stratified self-gravitating discs. 
The dispersion properties of modes in the presence of self- 
gravity as found here are one of the necessary things for 
studying mode propagation. 

Another point we want to raise concerns the spatial dis- 
tribution of temperature. In or der to realistically model the 
cooling of protoplanetary discs. [Bole v et akl l|2007h employed 
the radiative transfer technique. In the vertical direction, 
the radiative transfer equation was solved exactly assuming 
a plane-parallel atmosphere approximation, but in the ra- 
dial direction only the radiation diffusion approximation was 
employed. However, as our li near results (Fig. 11) and other 
non-linear simulations (e.g., iMeiia et al. I l2005l : iBolev et al.l 
I2OOQ) demonstrate, temperature and, therefore, opacity may 
vary on comparable scales in both the radial and vertical di- 
rections and have very complex structure in the non-linear 
regime. This implies that a more general radiative transfer 
treatment based on solving the ray equation in all directions. 
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rather than using the diffusion approximation in either di- 
rection, would be more appropriate for better understanding 
coohng processes. 
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